
<!DOCTYPE html>

<html>
  
<!-- Mirrored from docs.sympy.org/latest/modules/physics/mechanics/examples/lin_pend_nonmin_example.html by HTTrack Website Copier/3.x [XR&CO'2014], Sat, 15 Jan 2022 03:29:17 GMT -->
<!-- Added by HTTrack --><meta http-equiv="content-type" content="text/html;charset=utf-8" /><!-- /Added by HTTrack -->
<head>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="generator" content="Docutils 0.17.1: http://docutils.sourceforge.net/" />

    <title>Nonminimal Coordinates Pendulum &#8212; SymPy 1.9 documentation</title>
    <link rel="stylesheet" type="text/css" href="../../../../_static/pygments.css" />
    <link rel="stylesheet" type="text/css" href="../../../../_static/default.css" />
    <link rel="stylesheet" type="text/css" href="../../../../_static/graphviz.css" />
    <link rel="stylesheet" type="text/css" href="../../../../_static/plot_directive.css" />
    <link rel="stylesheet" type="text/css" href="../../../../../../live.sympy.org/static/live-core.css" />
    <link rel="stylesheet" type="text/css" href="../../../../../../live.sympy.org/static/live-autocomplete.css" />
    <link rel="stylesheet" type="text/css" href="../../../../../../live.sympy.org/static/live-sphinx.css" />
    
    <script data-url_root="../../../../" id="documentation_options" src="../../../../_static/documentation_options.js"></script>
    <script src="../../../../_static/jquery.js"></script>
    <script src="../../../../_static/underscore.js"></script>
    <script src="../../../../_static/doctools.js"></script>
    <script src="../../../../../../live.sympy.org/static/utilities.js"></script>
    <script src="../../../../../../live.sympy.org/static/external/classy.js"></script>
    <script src="../../../../../../live.sympy.org/static/live-core.js"></script>
    <script src="../../../../../../live.sympy.org/static/live-autocomplete.js"></script>
    <script src="../../../../../../live.sympy.org/static/live-sphinx.js"></script>
    <script async="async" src="../../../../../../cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest8331.js?config=TeX-AMS_HTML-full"></script>
    <script type="text/x-mathjax-config">MathJax.Hub.Config({"tex2jax": {"inlineMath": [["\\(", "\\)"]], "displayMath": [["\\[", "\\]"]]}})</script>
    
    <link rel="shortcut icon" href="../../../../_static/sympy-notailtext-favicon.ico"/>
    <link href="lin_pend_nonmin_example.html" rel="canonical" />
    
    <link rel="index" title="Index" href="../../../../genindex.html" />
    <link rel="search" title="Search" href="../../../../search.html" />
    <link rel="next" title="Examples for Physics/Mechanics" href="../examples.html" />
    <link rel="prev" title="Linearization in Physics/Mechanics" href="../linearize.html" /> 
  </head><body>
    <div class="related" role="navigation" aria-label="related navigation">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="../../../../genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="../../../../py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="../examples.html" title="Examples for Physics/Mechanics"
             accesskey="N">next</a> |</li>
        <li class="right" >
          <a href="../linearize.html" title="Linearization in Physics/Mechanics"
             accesskey="P">previous</a> |</li>
        <li class="nav-item nav-item-0"><a href="../../../../index.html">SymPy 1.9 documentation</a> &#187;</li>
          <li class="nav-item nav-item-1"><a href="../../../index.html" >SymPy Modules Reference</a> &#187;</li>
          <li class="nav-item nav-item-2"><a href="../../index.html" >Physics</a> &#187;</li>
          <li class="nav-item nav-item-3"><a href="../index.html" >Classical Mechanics</a> &#187;</li>
          <li class="nav-item nav-item-4"><a href="../linearize.html" accesskey="U">Linearization in Physics/Mechanics</a> &#187;</li>
        <li class="nav-item nav-item-this"><a href="#">Nonminimal Coordinates Pendulum</a></li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body" role="main">
            
  <section id="nonminimal-coordinates-pendulum">
<h1>Nonminimal Coordinates Pendulum<a class="headerlink" href="#nonminimal-coordinates-pendulum" title="Permalink to this headline">¶</a></h1>
<p>In this example we demonstrate the use of the functionality provided in
<a class="reference internal" href="../index.html#module-sympy.physics.mechanics" title="sympy.physics.mechanics"><code class="xref py py-mod docutils literal notranslate"><span class="pre">sympy.physics.mechanics</span></code></a> for deriving the equations of motion (EOM) for a pendulum
with a nonminimal set of coordinates. As the pendulum is a one degree of
freedom system, it can be described using one coordinate and one speed (the
pendulum angle, and the angular velocity respectively). Choosing instead to
describe the system using the <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> coordinates of the mass results in
a need for constraints. The system is shown below:</p>
<img alt="../../../../_images/pendulum_nonmin.svg" class="align-center" src="../../../../_images/pendulum_nonmin.svg" /><p>The system will be modeled using both Kane’s and Lagrange’s methods, and the
resulting EOM linearized. While this is a simple problem, it should illustrate
the use of the linearization methods in the presence of constraints.</p>
<section id="kane-s-method">
<h2>Kane’s Method<a class="headerlink" href="#kane-s-method" title="Permalink to this headline">¶</a></h2>
<p>First we need to create the <code class="docutils literal notranslate"><span class="pre">dynamicsymbols</span></code> needed to describe the system as
shown in the above diagram. In this case, the generalized coordinates <span class="math notranslate nohighlight">\(q_1\)</span> and
<span class="math notranslate nohighlight">\(q_2\)</span> represent the mass <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> coordinates in the inertial <span class="math notranslate nohighlight">\(N\)</span> frame.
Likewise, the generalized speeds <span class="math notranslate nohighlight">\(u_1\)</span> and <span class="math notranslate nohighlight">\(u_2\)</span> represent the velocities in
these directions. We also create some <code class="docutils literal notranslate"><span class="pre">symbols</span></code> to represent the length and
mass of the pendulum, as well as gravity and time.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">sympy.physics.mechanics</span> <span class="kn">import</span> <span class="o">*</span>
<span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">sympy</span> <span class="kn">import</span> <span class="n">symbols</span><span class="p">,</span> <span class="n">atan</span><span class="p">,</span> <span class="n">Matrix</span><span class="p">,</span> <span class="n">solve</span>
<span class="gp">&gt;&gt;&gt; </span><span class="c1"># Create generalized coordinates and speeds for this non-minimal realization</span>
<span class="gp">&gt;&gt;&gt; </span><span class="c1"># q1, q2 = N.x and N.y coordinates of pendulum</span>
<span class="gp">&gt;&gt;&gt; </span><span class="c1"># u1, u2 = N.x and N.y velocities of pendulum</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">q1</span><span class="p">,</span> <span class="n">q2</span> <span class="o">=</span> <span class="n">dynamicsymbols</span><span class="p">(</span><span class="s1">&#39;q1:3&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">q1d</span><span class="p">,</span> <span class="n">q2d</span> <span class="o">=</span> <span class="n">dynamicsymbols</span><span class="p">(</span><span class="s1">&#39;q1:3&#39;</span><span class="p">,</span> <span class="n">level</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">u1</span><span class="p">,</span> <span class="n">u2</span> <span class="o">=</span> <span class="n">dynamicsymbols</span><span class="p">(</span><span class="s1">&#39;u1:3&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">u1d</span><span class="p">,</span> <span class="n">u2d</span> <span class="o">=</span> <span class="n">dynamicsymbols</span><span class="p">(</span><span class="s1">&#39;u1:3&#39;</span><span class="p">,</span> <span class="n">level</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">L</span><span class="p">,</span> <span class="n">m</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">symbols</span><span class="p">(</span><span class="s1">&#39;L, m, g, t&#39;</span><span class="p">)</span>
</pre></div>
</div>
<p>Next, we create a world coordinate frame <span class="math notranslate nohighlight">\(N\)</span>, and its origin point <span class="math notranslate nohighlight">\(N^*\)</span>. The
velocity of the origin is set to 0. A second coordinate frame <span class="math notranslate nohighlight">\(A\)</span> is oriented
such that its x-axis is along the pendulum (as shown in the diagram above).</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Compose world frame</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">N</span> <span class="o">=</span> <span class="n">ReferenceFrame</span><span class="p">(</span><span class="s1">&#39;N&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pN</span> <span class="o">=</span> <span class="n">Point</span><span class="p">(</span><span class="s1">&#39;N*&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pN</span><span class="o">.</span><span class="n">set_vel</span><span class="p">(</span><span class="n">N</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>

<span class="gp">&gt;&gt;&gt; </span><span class="c1"># A.x is along the pendulum</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">theta1</span> <span class="o">=</span> <span class="n">atan</span><span class="p">(</span><span class="n">q2</span><span class="o">/</span><span class="n">q1</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span> <span class="o">=</span> <span class="n">N</span><span class="o">.</span><span class="n">orientnew</span><span class="p">(</span><span class="s1">&#39;A&#39;</span><span class="p">,</span> <span class="s1">&#39;axis&#39;</span><span class="p">,</span> <span class="p">[</span><span class="n">theta1</span><span class="p">,</span> <span class="n">N</span><span class="o">.</span><span class="n">z</span><span class="p">])</span>
</pre></div>
</div>
<p>Locating the pendulum mass is then as easy as specifying its location with in
terms of its x and y coordinates in the world frame. A <code class="docutils literal notranslate"><span class="pre">Particle</span></code> object is
then created to represent the mass at this location.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Locate the pendulum mass</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">P</span> <span class="o">=</span> <span class="n">pN</span><span class="o">.</span><span class="n">locatenew</span><span class="p">(</span><span class="s1">&#39;P1&#39;</span><span class="p">,</span> <span class="n">q1</span><span class="o">*</span><span class="n">N</span><span class="o">.</span><span class="n">x</span> <span class="o">+</span> <span class="n">q2</span><span class="o">*</span><span class="n">N</span><span class="o">.</span><span class="n">y</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pP</span> <span class="o">=</span> <span class="n">Particle</span><span class="p">(</span><span class="s1">&#39;pP&#39;</span><span class="p">,</span> <span class="n">P</span><span class="p">,</span> <span class="n">m</span><span class="p">)</span>
</pre></div>
</div>
<p>The kinematic differential equations (KDEs) relate the derivatives of the
generalized coordinates to the generalized speeds. In this case the speeds are
the derivatives, so these are simple. A dictionary is also created to map
<span class="math notranslate nohighlight">\(\dot{q}\)</span> to <span class="math notranslate nohighlight">\(u\)</span>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Calculate the kinematic differential equations</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">kde</span> <span class="o">=</span> <span class="n">Matrix</span><span class="p">([</span><span class="n">q1d</span> <span class="o">-</span> <span class="n">u1</span><span class="p">,</span>
<span class="gp">... </span>              <span class="n">q2d</span> <span class="o">-</span> <span class="n">u2</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">dq_dict</span> <span class="o">=</span> <span class="n">solve</span><span class="p">(</span><span class="n">kde</span><span class="p">,</span> <span class="p">[</span><span class="n">q1d</span><span class="p">,</span> <span class="n">q2d</span><span class="p">])</span>
</pre></div>
</div>
<p>The velocity of the mass is then the time derivative of the position from the
origin <span class="math notranslate nohighlight">\(N^*\)</span>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Set velocity of point P</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">P</span><span class="o">.</span><span class="n">set_vel</span><span class="p">(</span><span class="n">N</span><span class="p">,</span> <span class="n">P</span><span class="o">.</span><span class="n">pos_from</span><span class="p">(</span><span class="n">pN</span><span class="p">)</span><span class="o">.</span><span class="n">dt</span><span class="p">(</span><span class="n">N</span><span class="p">)</span><span class="o">.</span><span class="n">subs</span><span class="p">(</span><span class="n">dq_dict</span><span class="p">))</span>
</pre></div>
</div>
<p>As this system has more coordinates than degrees of freedom, constraints are
needed. The configuration constraints relate the coordinates to each other. In
this case the constraint is that the distance from the origin to the mass is
always the length <span class="math notranslate nohighlight">\(L\)</span> (the pendulum doesn’t get longer). Likewise, the velocity
constraint is that the mass velocity in the <code class="docutils literal notranslate"><span class="pre">A.x</span></code> direction is always 0 (no
radial velocity).</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="n">f_c</span> <span class="o">=</span> <span class="n">Matrix</span><span class="p">([</span><span class="n">P</span><span class="o">.</span><span class="n">pos_from</span><span class="p">(</span><span class="n">pN</span><span class="p">)</span><span class="o">.</span><span class="n">magnitude</span><span class="p">()</span> <span class="o">-</span> <span class="n">L</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">f_v</span> <span class="o">=</span> <span class="n">Matrix</span><span class="p">([</span><span class="n">P</span><span class="o">.</span><span class="n">vel</span><span class="p">(</span><span class="n">N</span><span class="p">)</span><span class="o">.</span><span class="n">express</span><span class="p">(</span><span class="n">A</span><span class="p">)</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">A</span><span class="o">.</span><span class="n">x</span><span class="p">)])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">f_v</span><span class="o">.</span><span class="n">simplify</span><span class="p">()</span>
</pre></div>
</div>
<p>The force on the system is just gravity, at point <code class="docutils literal notranslate"><span class="pre">P</span></code>.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Input the force resultant at P</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">R</span> <span class="o">=</span> <span class="n">m</span><span class="o">*</span><span class="n">g</span><span class="o">*</span><span class="n">N</span><span class="o">.</span><span class="n">x</span>
</pre></div>
</div>
<p>With the problem setup, the equations of motion can be generated using the
<code class="docutils literal notranslate"><span class="pre">KanesMethod</span></code> class. As there are constraints, dependent and independent
coordinates need to be provided to the class. In this case we’ll use <span class="math notranslate nohighlight">\(q_2\)</span> and
<span class="math notranslate nohighlight">\(u_2\)</span> as the independent coordinates and speeds:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Derive the equations of motion using the KanesMethod class.</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">KM</span> <span class="o">=</span> <span class="n">KanesMethod</span><span class="p">(</span><span class="n">N</span><span class="p">,</span> <span class="n">q_ind</span><span class="o">=</span><span class="p">[</span><span class="n">q2</span><span class="p">],</span> <span class="n">u_ind</span><span class="o">=</span><span class="p">[</span><span class="n">u2</span><span class="p">],</span> <span class="n">q_dependent</span><span class="o">=</span><span class="p">[</span><span class="n">q1</span><span class="p">],</span>
<span class="gp">... </span>                 <span class="n">u_dependent</span><span class="o">=</span><span class="p">[</span><span class="n">u1</span><span class="p">],</span> <span class="n">configuration_constraints</span><span class="o">=</span><span class="n">f_c</span><span class="p">,</span>
<span class="gp">... </span>                 <span class="n">velocity_constraints</span><span class="o">=</span><span class="n">f_v</span><span class="p">,</span> <span class="n">kd_eqs</span><span class="o">=</span><span class="n">kde</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="p">(</span><span class="n">fr</span><span class="p">,</span> <span class="n">frstar</span><span class="p">)</span> <span class="o">=</span> <span class="n">KM</span><span class="o">.</span><span class="n">kanes_equations</span><span class="p">([</span><span class="n">pP</span><span class="p">],[(</span><span class="n">P</span><span class="p">,</span> <span class="n">R</span><span class="p">)])</span>
</pre></div>
</div>
<p>For linearization, operating points can be specified on the call, or be
substituted in afterwards. In this case we’ll provide them in the call,
supplied in a list.  The <code class="docutils literal notranslate"><span class="pre">A_and_B=True</span></code> kwarg indicates to solve invert the
<span class="math notranslate nohighlight">\(M\)</span> matrix and solve for just the explicit linearized <span class="math notranslate nohighlight">\(A\)</span> and <span class="math notranslate nohighlight">\(B\)</span> matrices. The
<code class="docutils literal notranslate"><span class="pre">simplify=True</span></code> kwarg indicates to simplify inside the linearize call, and
return the presimplified matrices. The cost of doing this is small for simple
systems, but for larger systems this can be a costly operation, and should be
avoided.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Set the operating point to be straight down, and non-moving</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">q_op</span> <span class="o">=</span> <span class="p">{</span><span class="n">q1</span><span class="p">:</span> <span class="n">L</span><span class="p">,</span> <span class="n">q2</span><span class="p">:</span> <span class="mi">0</span><span class="p">}</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">u_op</span> <span class="o">=</span> <span class="p">{</span><span class="n">u1</span><span class="p">:</span> <span class="mi">0</span><span class="p">,</span> <span class="n">u2</span><span class="p">:</span> <span class="mi">0</span><span class="p">}</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">ud_op</span> <span class="o">=</span> <span class="p">{</span><span class="n">u1d</span><span class="p">:</span> <span class="mi">0</span><span class="p">,</span> <span class="n">u2d</span><span class="p">:</span> <span class="mi">0</span><span class="p">}</span>
<span class="gp">&gt;&gt;&gt; </span><span class="c1"># Perform the linearization</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span><span class="p">,</span> <span class="n">B</span><span class="p">,</span> <span class="n">inp_vec</span> <span class="o">=</span> <span class="n">KM</span><span class="o">.</span><span class="n">linearize</span><span class="p">(</span><span class="n">op_point</span><span class="o">=</span><span class="p">[</span><span class="n">q_op</span><span class="p">,</span> <span class="n">u_op</span><span class="p">,</span> <span class="n">ud_op</span><span class="p">],</span> <span class="n">A_and_B</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span>
<span class="gp">... </span>                             <span class="n">new_method</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> <span class="n">simplify</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span>
<span class="go">Matrix([</span>
<span class="go">[   0, 1],</span>
<span class="go">[-g/L, 0]])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">B</span>
<span class="go">Matrix(0, 0, [])</span>
</pre></div>
</div>
<p>The resulting <span class="math notranslate nohighlight">\(A\)</span> matrix has dimensions 2 x 2, while the number of total states
is <code class="docutils literal notranslate"><span class="pre">len(q)</span> <span class="pre">+</span> <span class="pre">len(u)</span> <span class="pre">=</span> <span class="pre">2</span> <span class="pre">+</span> <span class="pre">2</span> <span class="pre">=</span> <span class="pre">4</span></code>. This is because for constrained systems the
resulting <code class="docutils literal notranslate"><span class="pre">A_and_B</span></code> form has a partitioned state vector only containing
the independent coordinates and speeds. Written out mathematically, the system
linearized about this point would be written as:</p>
<div class="math notranslate nohighlight">
\[\begin{split}\begin{bmatrix} \dot{q_2} \\ \dot{u_2} \end{bmatrix} =
\begin{bmatrix} 0 &amp; 1 \\ \frac{-g}{L} &amp; 0 \end{bmatrix}
\begin{bmatrix} q_2 \\ u_2 \end{bmatrix}\end{split}\]</div>
</section>
<section id="lagrange-s-method">
<h2>Lagrange’s Method<a class="headerlink" href="#lagrange-s-method" title="Permalink to this headline">¶</a></h2>
<p>The derivation using Lagrange’s method is very similar to the approach using
Kane’s method described above. As before, we first create the
<code class="docutils literal notranslate"><span class="pre">dynamicsymbols</span></code> needed to describe the system. In this case, the generalized
coordinates <span class="math notranslate nohighlight">\(q_1\)</span> and <span class="math notranslate nohighlight">\(q_2\)</span> represent the mass <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> coordinates in the
inertial <span class="math notranslate nohighlight">\(N\)</span> frame.  This results in the time derivatives <span class="math notranslate nohighlight">\(\dot{q_1}\)</span> and
<span class="math notranslate nohighlight">\(\dot{q_2}\)</span> representing the velocities in these directions. We also create some
<code class="docutils literal notranslate"><span class="pre">symbols</span></code> to represent the length and mass of the pendulum, as well as
gravity and time.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">sympy.physics.mechanics</span> <span class="kn">import</span> <span class="o">*</span>
<span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">sympy</span> <span class="kn">import</span> <span class="n">symbols</span><span class="p">,</span> <span class="n">atan</span><span class="p">,</span> <span class="n">Matrix</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">q1</span><span class="p">,</span> <span class="n">q2</span> <span class="o">=</span> <span class="n">dynamicsymbols</span><span class="p">(</span><span class="s1">&#39;q1:3&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">q1d</span><span class="p">,</span> <span class="n">q2d</span> <span class="o">=</span> <span class="n">dynamicsymbols</span><span class="p">(</span><span class="s1">&#39;q1:3&#39;</span><span class="p">,</span> <span class="n">level</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">L</span><span class="p">,</span> <span class="n">m</span><span class="p">,</span> <span class="n">g</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">symbols</span><span class="p">(</span><span class="s1">&#39;L, m, g, t&#39;</span><span class="p">)</span>
</pre></div>
</div>
<p>Next, we create a world coordinate frame <span class="math notranslate nohighlight">\(N\)</span>, and its origin point <span class="math notranslate nohighlight">\(N^*\)</span>. The
velocity of the origin is set to 0. A second coordinate frame <span class="math notranslate nohighlight">\(A\)</span> is oriented
such that its x-axis is along the pendulum (as shown in the diagram above).</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Compose World Frame</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">N</span> <span class="o">=</span> <span class="n">ReferenceFrame</span><span class="p">(</span><span class="s1">&#39;N&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pN</span> <span class="o">=</span> <span class="n">Point</span><span class="p">(</span><span class="s1">&#39;N*&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pN</span><span class="o">.</span><span class="n">set_vel</span><span class="p">(</span><span class="n">N</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="c1"># A.x is along the pendulum</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">theta1</span> <span class="o">=</span> <span class="n">atan</span><span class="p">(</span><span class="n">q2</span><span class="o">/</span><span class="n">q1</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span> <span class="o">=</span> <span class="n">N</span><span class="o">.</span><span class="n">orientnew</span><span class="p">(</span><span class="s1">&#39;A&#39;</span><span class="p">,</span> <span class="s1">&#39;axis&#39;</span><span class="p">,</span> <span class="p">[</span><span class="n">theta1</span><span class="p">,</span> <span class="n">N</span><span class="o">.</span><span class="n">z</span><span class="p">])</span>
</pre></div>
</div>
<p>Locating the pendulum mass is then as easy as specifying its location with in
terms of its x and y coordinates in the world frame. A <code class="docutils literal notranslate"><span class="pre">Particle</span></code> object is
then created to represent the mass at this location.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Create point P, the pendulum mass</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">P</span> <span class="o">=</span> <span class="n">pN</span><span class="o">.</span><span class="n">locatenew</span><span class="p">(</span><span class="s1">&#39;P1&#39;</span><span class="p">,</span> <span class="n">q1</span><span class="o">*</span><span class="n">N</span><span class="o">.</span><span class="n">x</span> <span class="o">+</span> <span class="n">q2</span><span class="o">*</span><span class="n">N</span><span class="o">.</span><span class="n">y</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">P</span><span class="o">.</span><span class="n">set_vel</span><span class="p">(</span><span class="n">N</span><span class="p">,</span> <span class="n">P</span><span class="o">.</span><span class="n">pos_from</span><span class="p">(</span><span class="n">pN</span><span class="p">)</span><span class="o">.</span><span class="n">dt</span><span class="p">(</span><span class="n">N</span><span class="p">))</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pP</span> <span class="o">=</span> <span class="n">Particle</span><span class="p">(</span><span class="s1">&#39;pP&#39;</span><span class="p">,</span> <span class="n">P</span><span class="p">,</span> <span class="n">m</span><span class="p">)</span>
</pre></div>
</div>
<p>As this system has more coordinates than degrees of freedom, constraints are
needed. In this case only a single holonomic constraints is needed: the
distance from the origin to the mass is always the length <span class="math notranslate nohighlight">\(L\)</span> (the pendulum
doesn’t get longer).</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Holonomic Constraint Equations</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">f_c</span> <span class="o">=</span> <span class="n">Matrix</span><span class="p">([</span><span class="n">q1</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">q2</span><span class="o">**</span><span class="mi">2</span> <span class="o">-</span> <span class="n">L</span><span class="o">**</span><span class="mi">2</span><span class="p">])</span>
</pre></div>
</div>
<p>The force on the system is just gravity, at point <code class="docutils literal notranslate"><span class="pre">P</span></code>.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Input the force resultant at P</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">R</span> <span class="o">=</span> <span class="n">m</span><span class="o">*</span><span class="n">g</span><span class="o">*</span><span class="n">N</span><span class="o">.</span><span class="n">x</span>
</pre></div>
</div>
<p>With the problem setup, the Lagrangian can be calculated, and the equations of
motion formed. Note that the call to <code class="docutils literal notranslate"><span class="pre">LagrangesMethod</span></code> includes the
Lagrangian, the generalized coordinates, the constraints (specified by
<code class="docutils literal notranslate"><span class="pre">hol_coneqs</span></code> or <code class="docutils literal notranslate"><span class="pre">nonhol_coneqs</span></code>), the list of (body, force) pairs, and the
inertial frame. In contrast to the <code class="docutils literal notranslate"><span class="pre">KanesMethod</span></code> initializer, independent and
dependent coordinates are not partitioned inside the <code class="docutils literal notranslate"><span class="pre">LagrangesMethod</span></code>
object. Such a partition is supplied later.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Calculate the lagrangian, and form the equations of motion</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">Lag</span> <span class="o">=</span> <span class="n">Lagrangian</span><span class="p">(</span><span class="n">N</span><span class="p">,</span> <span class="n">pP</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">LM</span> <span class="o">=</span> <span class="n">LagrangesMethod</span><span class="p">(</span><span class="n">Lag</span><span class="p">,</span> <span class="p">[</span><span class="n">q1</span><span class="p">,</span> <span class="n">q2</span><span class="p">],</span> <span class="n">hol_coneqs</span><span class="o">=</span><span class="n">f_c</span><span class="p">,</span> <span class="n">forcelist</span><span class="o">=</span><span class="p">[(</span><span class="n">P</span><span class="p">,</span> <span class="n">R</span><span class="p">)],</span> <span class="n">frame</span><span class="o">=</span><span class="n">N</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">lag_eqs</span> <span class="o">=</span> <span class="n">LM</span><span class="o">.</span><span class="n">form_lagranges_equations</span><span class="p">()</span>
</pre></div>
</div>
<p>Next, we compose the operating point dictionary, set in the hanging at rest
position:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Compose operating point</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">op_point</span> <span class="o">=</span> <span class="p">{</span><span class="n">q1</span><span class="p">:</span> <span class="n">L</span><span class="p">,</span> <span class="n">q2</span><span class="p">:</span> <span class="mi">0</span><span class="p">,</span> <span class="n">q1d</span><span class="p">:</span> <span class="mi">0</span><span class="p">,</span> <span class="n">q2d</span><span class="p">:</span> <span class="mi">0</span><span class="p">,</span> <span class="n">q1d</span><span class="o">.</span><span class="n">diff</span><span class="p">(</span><span class="n">t</span><span class="p">):</span> <span class="mi">0</span><span class="p">,</span> <span class="n">q2d</span><span class="o">.</span><span class="n">diff</span><span class="p">(</span><span class="n">t</span><span class="p">):</span> <span class="mi">0</span><span class="p">}</span>
</pre></div>
</div>
<p>As there are constraints in the formulation, there will be corresponding
Lagrange Multipliers. These may appear inside the linearized form as well, and
thus should also be included inside the operating point dictionary.
Fortunately, the <code class="docutils literal notranslate"><span class="pre">LagrangesMethod</span></code> class provides an easy way of solving
for the multipliers at a given operating point using the <code class="docutils literal notranslate"><span class="pre">solve_multipliers</span></code>
method.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="c1"># Solve for multiplier operating point</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">lam_op</span> <span class="o">=</span> <span class="n">LM</span><span class="o">.</span><span class="n">solve_multipliers</span><span class="p">(</span><span class="n">op_point</span><span class="o">=</span><span class="n">op_point</span><span class="p">)</span>
</pre></div>
</div>
<p>With this solution, linearization can be completed. Note that in contrast to
the <code class="docutils literal notranslate"><span class="pre">KanesMethod</span></code> approach, the <code class="docutils literal notranslate"><span class="pre">LagrangesMethod.linearize</span></code> method also
requires the partitioning of the generalized coordinates and their time
derivatives into independent and dependent vectors.  This is the same as what
was passed into the <code class="docutils literal notranslate"><span class="pre">KanesMethod</span></code> constructor above:</p>
<div class="doctest highlight-default notranslate"><div class="highlight"><pre><span></span><span class="gp">&gt;&gt;&gt; </span><span class="n">op_point</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">lam_op</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="c1"># Perform the Linearization</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span><span class="p">,</span> <span class="n">B</span><span class="p">,</span> <span class="n">inp_vec</span> <span class="o">=</span> <span class="n">LM</span><span class="o">.</span><span class="n">linearize</span><span class="p">([</span><span class="n">q2</span><span class="p">],</span> <span class="p">[</span><span class="n">q2d</span><span class="p">],</span> <span class="p">[</span><span class="n">q1</span><span class="p">],</span> <span class="p">[</span><span class="n">q1d</span><span class="p">],</span>
<span class="gp">... </span>                            <span class="n">op_point</span><span class="o">=</span><span class="n">op_point</span><span class="p">,</span> <span class="n">A_and_B</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">A</span>
<span class="go">Matrix([</span>
<span class="go">[     0, 1],</span>
<span class="go">[-g/L, 0]])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">B</span>
<span class="go">Matrix(0, 0, [])</span>
</pre></div>
</div>
<p>The resulting <span class="math notranslate nohighlight">\(A\)</span> matrix has dimensions 2 x 2, while the number of total states
is <code class="docutils literal notranslate"><span class="pre">2*len(q)</span> <span class="pre">=</span> <span class="pre">4</span></code>. This is because for constrained systems the resulting
<code class="docutils literal notranslate"><span class="pre">A_and_B</span></code> form has a partitioned state vector only containing the independent
coordinates and their derivatives. Written out mathematically, the system
linearized about this point would be written as:</p>
<div class="math notranslate nohighlight">
\[\begin{split}\begin{bmatrix} \dot{q_2} \\ \ddot{q_2} \end{bmatrix} =
\begin{bmatrix} 0 &amp; 1 \\ \frac{-g}{L} &amp; 0 \end{bmatrix}
\begin{bmatrix} q_2 \\ \dot{q_2} \end{bmatrix}\end{split}\]</div>
</section>
</section>


            <div class="clearer"></div>
          </div>
        </div>
      </div>
      <div class="sphinxsidebar" role="navigation" aria-label="main navigation">
        <div class="sphinxsidebarwrapper">
            <p class="logo"><a href="../../../../index.html">
              <img class="logo" src="../../../../_static/sympylogo.png" alt="Logo"/>
            </a></p>
  <h3><a href="../../../../index.html">Table of Contents</a></h3>
  <ul>
<li><a class="reference internal" href="#">Nonminimal Coordinates Pendulum</a><ul>
<li><a class="reference internal" href="#kane-s-method">Kane’s Method</a></li>
<li><a class="reference internal" href="#lagrange-s-method">Lagrange’s Method</a></li>
</ul>
</li>
</ul>

  <h4>Previous topic</h4>
  <p class="topless"><a href="../linearize.html"
                        title="previous chapter">Linearization in Physics/Mechanics</a></p>
  <h4>Next topic</h4>
  <p class="topless"><a href="../examples.html"
                        title="next chapter">Examples for Physics/Mechanics</a></p>
  <div role="note" aria-label="source link">
    <h3>This Page</h3>
    <ul class="this-page-menu">
      <li><a href="../../../../_sources/modules/physics/mechanics/examples/lin_pend_nonmin_example.rst.txt"
            rel="nofollow">Show Source</a></li>
    </ul>
   </div>
<div id="searchbox" style="display: none" role="search">
  <h3 id="searchlabel">Quick search</h3>
    <div class="searchformwrapper">
    <form class="search" action="https://docs.sympy.org/latest/search.html" method="get">
      <input type="text" name="q" aria-labelledby="searchlabel" autocomplete="off" autocorrect="off" autocapitalize="off" spellcheck="false"/>
      <input type="submit" value="Go" />
    </form>
    </div>
</div>
<script>$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related" role="navigation" aria-label="related navigation">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="../../../../genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="../../../../py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="../examples.html" title="Examples for Physics/Mechanics"
             >next</a> |</li>
        <li class="right" >
          <a href="../linearize.html" title="Linearization in Physics/Mechanics"
             >previous</a> |</li>
        <li class="nav-item nav-item-0"><a href="../../../../index.html">SymPy 1.9 documentation</a> &#187;</li>
          <li class="nav-item nav-item-1"><a href="../../../index.html" >SymPy Modules Reference</a> &#187;</li>
          <li class="nav-item nav-item-2"><a href="../../index.html" >Physics</a> &#187;</li>
          <li class="nav-item nav-item-3"><a href="../index.html" >Classical Mechanics</a> &#187;</li>
          <li class="nav-item nav-item-4"><a href="../linearize.html" >Linearization in Physics/Mechanics</a> &#187;</li>
        <li class="nav-item nav-item-this"><a href="#">Nonminimal Coordinates Pendulum</a></li> 
      </ul>
    </div>
    <div class="footer" role="contentinfo">
        &#169; Copyright 2021 SymPy Development Team.
      Last updated on Sep 30, 2021.
      Created using <a href="https://www.sphinx-doc.org/">Sphinx</a> 4.1.2.
    </div>
  </body>

<!-- Mirrored from docs.sympy.org/latest/modules/physics/mechanics/examples/lin_pend_nonmin_example.html by HTTrack Website Copier/3.x [XR&CO'2014], Sat, 15 Jan 2022 03:29:18 GMT -->
</html>